Numerical and analytical investigation of the free 
boundary confluence for the phase field system. 
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Abstract 

In this paper we numerically research the solutions of the phase field system for 
the spherically symmetric Stefan-Gibbs-Thomson problem in the case of interaction 



1 INTRODUCTION. STATEMENT OF THE PROB- 



The main goal of this paper is the numerical research of the effect of the soliton type distur- 
t— ( ■ bance of the temperature in the point of the contact of the free boundaries. This effect we 
consider in the case of the phase field model for the Stefan-Gibbs-Thomson problem. The 
', difference between this problem and the classical Stefan problem is that the surface tension 
is taken into account in the Stefan-Gibbs-Thomson problem. At the beginning of the paper 
we briefly give the main analytic results about the confluence of the free boundaries and 
^ ■ then we illustrate these results by computer simulation. 

The effect of the soliton type (negative) disturbance of the temperature a of the free 
T"! \ boundaries in the point of the contact is shown in Fig. [TJ 

. ^ ! We present the results of the computer simulation for the phase field system in the 

^| spherically symmetric case. Namely, we consider the domain Q = f2 x [0,£i], where f2 = 

[Ri, R%\ is the spherical layer in the spherical coordinates. We assume that the domain Q is 

divided into the three layers Qf 2 (t) an d as follows: 

Qf(t) = [R 1 ,r 1 (t)], fi-(t) = [n(t),r a (t)], Q+(t) = [r 2 (t), R 2 ], 

where r»(t) = % — 1, 2 is the free boundaries between phases "+" and "— ". We assume 

that the phase "+" occupies the layer flf 2 (t), and the phase "— " occupies the layer Q~(t) 
(see Fig. EJ. 

In this case the phase field system have the form pQ 

du u — u 3 

L9 = - — , eLu x9 = 0, 1 

dt e w 
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Figure 1: Effect of the soliton type disturbance of the temperature (e = 0.003). a) Coordinate 
and time dependence of the temperature, b) Soliton type disturbance of the temperature in 
the neighborhood of the point of the contact of the free boundaries, c) Initial data of the 
temperature and the evolution of the free boundaries ri(t), ^(t) in prime. 



where x = y/2/3, 



The function 6 = 6(r,t,e) has the meaning of the temperature, and the function u = 
u(r,t,e) (so-called the order function) determines the phase state of the medium. Namely, 
the value u ~ — 1 corresponds to the phase "— " in the layer Q~(t), and the value u ~ 1 
corresponds to the the phase "+" in the layers f^ 2 (£)- 

We denote 

M-^ (2) 
The phase field system ([T]) we can write in the form 

da d 2 a du 

~dt~d^ = ~ r dt' ( ' 

u - u 3 a 

eLu = h x-. (4) 

e r 

Passing to the limit as e — > in PJ), Ji) we obtain the Stefan-Gibbs-Thomson problem 




Figure 2: The free boundaries in the spherically symmetric case. 



for the each free boundary (see [U El E] ) 
da d 2 a 



re[Rx,R 2 }, r^fi(t), i 
2 



dt dr 2 
da 



1,2, 



dr 



1 



2r i (t)f' i (t). 



(5) 
(6) 

(7) 



r=fi(t) 



This passage to the limit is possible, for example, in the case where the corresponding limit 
problems have classical solutions. In this case, the weak limits as e — > of solutions ([3]), (j3J) 
give these solutions [H El E] ■ The existence of the classical solution a of problem (j5])-(I7j) is 
discussed in [8]. 

By t* G (0, ti) we denote the instant of time of the confluence of the free boundaries, and 
r* = ri(t*) = f 2 (t*) is the sphere of the contact of the free boundaries. 

The smooth approximations (i.e. the approximate solution of the system P| with 
misalignment that is small in the weak sense as e — > 0, see [2]) of the Stefan- Gibbs-Thomson 
problem (of the phase field system) is constructed in [3J in the one dimensional case and this 
smooth approximations exists in the time included the instant of the contact t* . Moreover, 
this smooth approximations is constructed on the assumption of the existence of the classical 
solution of the limit problem as t ^ t* — 5, where 5 > is any number. This approximations 
admits the passage to the limit as e — > 0. 

It is possible to show [2] that in the common sense the asymptotic solution of the system 
(El), (H has the form 



a-{r,t) + (a + (r,t) - cr"(r,t)) 



f 2 (t) - r 



r — fi(t) 



u 



fi(t) — r 



(8) 
(9) 



V — ro(t) 



r — ri(t) r — f^if) 



as t ^ t* - 5, 6 > 0. Here u^z) ^0,lasz^ ^oo, u[ k \z) G S(R*) as fc > 0, r,(t), z = 1, 2 
is the smooth functions, r\ ^ r 2 , Wo(<£) — tanh(z), and u(t,z 1 ,z 2 ) G C°°([0,t*]; S(R^)). By 
S(]R n ) we denote the Schwartz space of smooth rapidly decreases functions. If the initial 
data for ([3]), (j4j) has the form (E]), (EJ) at t — 0, then, for £ ^ t* — 5, we have the estimate 



\u — u c 



C(0,T; L 2 



|a-< s ;£ 2 (Q)|| < 



/i ^ 3/2, 



where (a, u) is a solution of system (E]),(j4j) (see [2] and references in). Here Q = ft x [0, t* — S), 
and the constant c is independent of e. 

The main obstacle to the construction of approximations of solution in the case of con- 
fluence of free boundaries is the fact that, instead of an ordinary differential equation whose 
solution is the function uq(z) |U [2], in the case of confluence of free boundaries, we must 
deal with a partial differential equation for which the explicit form of the exact solution is 
unknown. 

In papers [31 0] using the the weak asymptotic method the solution of the phase field 
system was constructed that describes the confluence of the free boundaries. Namely, the 
ansatz of the order function has the form. 



^*2 ^2 



(10) 



wo 



^•2 ^2 



w /3 



r — rt 



The temperature is sought in the form 

a = e(r)T + q, 



e(r) G Cq°([Ri, R 2 )), e = 1 for r G [fi(0), f 2 (0)]. Here T is the model of the temperature, i.e. 
it is the function of the simplest structure which describes the behaviors of the temperature 
qualitatively correct, and q is an unknown smooth function. Namely, 



Here 
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f = 7i"(t)(ri - r)H(ri - r) + 7^ (i)(r - r 2 )H(r - r 2 ) 

r 2 - ri 
r 2 - n 



= 7i-7 2 (2 ^ 7i+7 2 
2 1 

7i + 72 



(12) 



2 

fcl - &2 



-(r 2 



r 



— (r — r 



; >n 7i ~ 72 

,2\ fc l + ^2 



fei = rir u + 2, fc 2 = r 2 r 2t + 2, 



^(t,e) = r 2 (t, £ ) -rl(t,e), 



(13) 
(14) 



So, from the given above formulas we see that the model T is linear in r in the layers Qf(t), 
% = 1, 2 and parabolic in r in the layer Q~(t). 

In paper [3] the formulas are obtained those determine the functions contained in ansatzes 
(TTO]) . (TTTT) in the one dimensional case. 
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Figure 3: Confluence of the free boundaries in the case of the curvatures of the free boundaries 
have the unlike signs. 

The analysis of these formulas shows that the contained in the phase field system tem- 
perature 9 has the soliton type disturbance in the neighborhood of the point of the contact. 
This disturbance is localized in the space coordinates and in time. The "width" of this lo- 
calization is proportional to e. In [3] the limit as e — > of the amplitude of the disturbance 
is derived. In the one dimensional case the value of this amplitude is given by formula 

\9]\ ^ _ = — lim lim — !— — - , 

l l\x=x*, t=t* t ^t*-0e^0 2 

where V\ and V 2 are the velocities of the merging (one-dimensional) free boundaries. 

The analytic treatment based on the weak asymptotics method shows that such effect is 
general. For example, the exactly same formula is correct in the case of merging globe layers 
[4] those is shown in Figj2j Namely, 

\9\ I . ^ = — lim hm lw — — , 

L J \r=r*, t=t* 2 ' 

V\ n and V 2n are the normal velocities of the merging three-dimensional free boundaries (cor- 
respondingly, Ti(t) and T 2 (t)). If the free boundaries are asymmetrical, then their principle 
curvatures have the like signs at the point of the contact. In this case the amplitude of the 
jump of the temperature is determined by formula 



( 




+ 


v 2n 


V 2 



\0}\, w .,. w =-limlim 1 ln| "' 2nl -\Jd-JC 

l J I i,»,2 = i',«*r , t=t* +^/*_n c -^n V 1 



\{x,y,z)=(x*,y*,z*), t=t* t ^ t *_o £ ^. 

in the instant of the confluence and at the point of the contact. Here 9 = 9(x,y, z,t), 
(x*,y*,z*) is the point of the contact of the free boundaries, and /Ci, /C 2 are the principle 
curvatures of the free boundaries at the point of the contact. 

Let us consider another situation. We assume that the solution of our problem is sym- 
metric about y-axe (see FigH]). By yi(t) and y 2 (t) we denote the points in which the free 
boundaries cross the y-axe and these points is situated on the shorts distance one from an- 
other. We assume that M* = (0, y*, 0) is the point of the contact of the free boundaries. The 
principle difference between this situation and the situation in Figj2] is that in the instant 
of the contact t* the principle curvatures of the free boundaries have the unlike signs at 
the point of the contact M*. In this case the amplitude of the jump of the temperature is 
determined by formula 

wn .nw *. = - lim lim ( 1W 2 " + |/C 5 

L J \(x,y,z)=(p,y*,0), t=t* t -rt.-o £ ^0\ 2 

The mentioned soliton type disturbance is observed in the numerical experiments, see 
FigHEOj 



We note that beside the our own papers we know only the single paper of A.M. Meirmanov 
and B.A. Zaltsman [TJ. In this paper the problem of the confluence of the free boundaries 
in the case of the Hele-Shaw problem is derived and this problem is the special case of the 
Stefan- Gibbs-Thomson problem. The regularized (not limit) problem is considered in [7J 
and in this case the effect of the soliton type disturbance is not shown. The analog of this 
fact is the following problem. Let us consider the heat equation in a rectilinear segment and 
with not zero (= 1) Dirichlet's initial condition. Clearly, the solution of this problem is zero. 
However, the numerical solution is not identical zero for the different scheme with a node in 
the point x — 0. 



2 DIFFERENT SCHEME 

The choice of the different scheme for system (J3j) , (SJ) is based on some ideas that are 
sufficiently general for solving nonlinear equations. Namely, we calculate the heat equation 
([3]) at the next time step and, consequently split system 0, (@J. We associate Eq. OH) with 
implicit different scheme 

n n k+1 

-r-s£<S- 2£ £^= as) 

1 i ' i 

i=l i=l 

= £ _ (( w *)3 + 3 ( w *)2( w *+l _ „*))] + — • 

i=l 

Here u k is the mesh order function at the fc-th time step, a k is the mesh "temperature" at 
the fc-th time step. 

The first term (the term in the square brackets) in the right hand part of Eq. (fT5l) is 
obtained as following. We denote 

F{u) = u - u 3 . 

We linearize the mesh function F{u k+l ) as following 

F(u k+1 ) w F(u k ) + F' u (u k )(u k+l - u k ). 

System (fT5l) is completed with initial and boundary conditions in the points r = Rj, 
j = 1,2. It is clear that Eqs. f|T5l) is the three-point equations relatively to u k+l and these 
equations are solved by the sweep method. 

To calculate the mesh "temperature" a k+1 we use the standard different scheme for the 
heat equation Q. Namely, 



n 

8=1 1=1 



r-E^ = -E^r- as) 



Equations ffTB]) is also solved by the sweep method. 

We use the main segment r e [1, 2] (i.e., Ri = 1, R<i = 2) for the numerical simulation of 
the process of the confluence. 

The initial data we choose as provided by the structure of the weak asymptotic solution 
Q). We use the function 

u° = 1 + tanh ^ + tanh 2 - (17) 




Figure 4: The initial data for the order function u for the different values of e: 1. e = 0.025, 
2. e = 0.01, 3. s = 0.007, 4. s = 0.005, 5. s = 0.003. 




Figure 5: The initial data for the temperature a. The initial positions of the free boundaries 
are r° = 1.25, r° 2 = 1.75 

as the initial data for the order function u. Here r° = 1.25, r° = 1.75 determine the initial 
position of the free boundaries (see Fig j4l5l . 

The initial data a for the "temperature" is taken in the form in Figj5] as provided by 
model (fl2"j) . Namely, the function o° is parabolic in the domain with phase "— " (between 
the free boundaries), and the function cr° is linear in the domain with phase "+". At the 
same time, the dependence between the initial positions and the initial velocities of the free 
boundaries and the values of a is determined by formula ([6]) in the boundary points r°, r°- 

It is clear that the considered initial data differ from the exact solution of the problem. 
Nevertheless it is known that the solutions of problem ([3]), (j4j) converge to solutions of the 
limit Stefan-Gibbs-Thomson problem as e — > and on the sufficiently general assumptions. 
More other, the solution of the Cauchi problem come to self-similar regime by the widely 
known behaviours of the semi-linear parabolic equations (see [5]). Beside it, the formulas for 
the asymptotic solution are the deformations of the formulas for the semi-similar solution. 
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tmin 




Gymim "mini & ) 


0.025 


0.08432 


1 1 1 o 

1.448 


1.373699 


0.01 


0.09690 


1.427 


1.051388 


0.007 


0.10016 


1.422 


0.967825 


0.005 


0.10260 


1.419 


0.902497 


0.003 


0.10549 


1.415 


0.587919 



Table 1: Coordinate and time dependence of the jump of the temperature a. r min and t 
are values of the coordinate r and time t at which the jump becomes minimum. 



3 THE RESULTS OF NUMERICAL SIMULATION 

Let us consider the results of the numerical simulation of the process of the confluence of 
the free boundaries. 

Some graphics given below illustrate the common behaviors of the solution of system 
([3]), (j3J) and confirm the certainty of the numerical results. We simulate for the different 
(decreasing) values of the parameter e, e — 0.025, e = 0.01, e = 0.007, 4. e — 0.005, 
e = 0.003. The mesh is equal h = 10~ 3 and the time step is equal r = 10~ 5 . 

In FigJHl El are shown the profiles of the temperature a and the order function u in the 
neighborhood of the point of the contact of the free boundaries (for fixed r* pa r = 1.42) for 
different values of the parameter e. 

From the graphics in Figj6]we see that the "width" of the neighborhood of the instant of 
time of the jump of the temperature decrease (it is proportional to e) owing to decreasing e. 
At the same time, the jump occurs at the the different instant of time for the different values 
of the parameter e. Namely, in dependence of the decreasing of e the sequence of the instant 
of time tmi n (i n which the jump becomes minimum) increases, see also Table [TJ Beside it, 
from Table [1] we see that the intervals grow short between the neighboring instants t min as 
the parameter e decreases. The cause of this fact is that the width of the transition zone 
of the order function u is different for different values of e, see formula (IT7|) . Namely, the 
width of transition zone is proportional e (as larger the value of e, as wider the transition 
zone), see Fig. [Hand [TJ [2]. As a result of this fact the contact of the transition zones (and, 
consequently the beginning of the confluence of the free boundaries) occurs previously for 
the simulation process with larger e, see Fig. |SJ Clearly, the width of the confluence is also 
proportional to e respect to r. From Fig. [S]we see that for the fixed instant of time t = 0.082 
the free boundaries corresponding to the graphic 1 (e = 0.025) complete the confluence, but 
at the same time the free boundaries corresponding another graphics (for smaller e) are on 
the sufficiently large distance. 

From Fig. [7] we see that as smaller e as quickly transit occurs from the phase "— " to the 
phase "+". It is clear that u(r*,t) —> sign(t — t min ) as e — > 0. 

Beside it, from Figj6]we see that the minimum of the jump for the graphic 5 (e = 0.003) 
is smaller then the minimum of the jump for the graphic 4 (e = 0.005). But in real this is 
not true and we observe reversed dependence, see Table [1] and Figj9j Here are shown the 
graphics of the temperature a in the fixed point r = r min , where r min is the value of the 
coordinate r, in which the jump of the temperature becomes zero. This fact means that the 
shift of the minimum of the jump of the temperature a depend on e in the coordinate r, see 
Table [TJ From Table [TJ we see that this shift is not lager the difference between the width 
of the transition zones (difference between the the corresponding values of the parameter e) 
for the corresponding graphics, see Figj9) 
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Figure 6: The jump of the temperature a in the neighborhood of the point of the contact of 
the free boundaries (r* « r = 1.42): 1. e = 0.025, 2. e = 0.01, 3. e = 0.007, 4. e = 0.005, 5. 
£ = 0.003. 
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Figure 7: The order function u in the neighborhood of the point of the contact of the free 
boundaries (r* w 1.42): 1. e = 0.025, 2. e = 0.01, 3. e = 0.007, 4. e = 0.005, 5. e = 0.003. 





Figure 9: The jump of the temperature a: 1. e = 0.025, r = 1.448, 2. e = 0.01, r = 1.427, 
3. e = 0.007, r = 1.422, 4. e = 0.005, r = 1.419, 5. £ = 0.003, r = 1.415. 



in 




Figure 10: The dynamic of the jump of the temperature a for e = 0.003: 1. t = 0.1052, 2. 
t = 0.10539, 3. t = 0.10544, 4. t = t min = 0.10549. 

From Tabled] we see that the intervals between the neighborhood values of i m j n decrease 
as e decreases. More other, this intervals is the order of e. At the same time the distance 
between the neighborhood points r min grows short as e decreases. 

The dynamic of the jump of the temperature at time is shown in Fig. [10] for e = 0.003. 

So we can conclude that the process of the confluence of the free boundaries is local- 
ized in the coordinate and at the time. Consequently, the soliton type disturbance of the 
temperature is localized in the coordinate and at the time (the width of the localization is 
proportional toe), see Fig. CD 

We denote that the values of the parameter e > 0.003 is not possible in the our simulation. 
For e = 0.003 and h = 0.001 we obtain that the three nodes (the minimal number of the 
nodes that necessary to calculate the second difference derivation, see) of the mesh appear 
in the transition zone, see (1151) . (TIB]) . The computer simulation shows that the numerical 
solution is unstable in the neighborhood of the confluence of the free boundaries for e = 0.002 
and the different scheme do not give solution for e = 0.001. 

The series of the graphics demonstrates the stability of the numerical solution of system 
(j3J), (jl]) respect to e. In FigJTT] [12] are shown the graphics of the dependence of the order 
function u and of the temperature a on time in the fixed point r = 1.15. We note that the 
analogous results are correct for the another value of r except the point of the contact of the 
free boundary. 

From the graphics in Fig{TT]we see that at the beginning (t ^ 0.01) the solution u undergo 
the sharp jump and thereafter stables (undergo to near semi-similar regime). At the same 
time, this stabilization is wavelike. More other, the initial interval of the time in which the 
solution undergo the most strong change (jump) grows short as e decreases. More other the 
distance between the neighborhood graphics grows short as e decreases and this distance is 
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Figure 11: The order function u in r = 1.15: I. e = 0.025, 2. e = 0.01, 3. e = 0.007, 4. 
£ = 0.005, 5. e = 0.003. 
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Figure 12: The temperature a in r = 1.15: 1. e = 0.025, 2. e = 0.01, 3. e = 0.007, 4. 
e = 0.005, 5. e = 0.003. 
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Figure 14: The dynamic of the free boundary r 2 (t) for the different values of the parameter 
e: 1. e = 0.025, 2. e = 0.01, 3. e = 0.007, 4. e = 0.005, 5. e = 0.003. 

proportional. 

From the graphics in Figjl2] we see that the temperature a is stable respect to the 
variation of e. For t > 0.09 the difference between the graphics due to the process of the 
confluence for the graphics for the larger e begins early than for the graphics for the smaller e. 

The last series of the numerical experiments deals with the verification of the analytical 
formula for the amplitude of the jump of the temperature. Namely, in paper j3] is obtained 
the formula for the amplitude of the jump 



or taken into account (T5]) we have 



'10 '20 

t=t* _ O ' 



MU. = - Hm lim ^MzMiMM , (18 ) 



In FigJTBl [TUthe dependence of the free boundaries ri(t) and ^(t) on time is shown for 
the different values of the parameter e. We see that the functions rx(t) and ^(t) is linear in 
the initial interval of the shifting of the free boundaries. The nonlinearity of the confluence 
depends on the shifting of the free boundary near the point of the contact and the lines are 
distorted. We see that as smaller e as smaller the neighborhood of the instant t* in which 
the distortion is sensed. 
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ri 


r[ 


r 2 


/ 


r 1 
[&an\ 


r i 
[CTcal\ 


0.025 


1 on a 

1.394 


1.4 


1 A O f 

1.485 


-3.2 


o n 1 A A o 

3.01448 


1.5 


n n i 

0.01 


i /inn 

1.402 


1 


1 a a n 

1.446 


-3.2 


3.3518 


1.7 


0.007 


1.403 


0.8 


1.438 


-3.2 


2.86 


1.9 


0.005 


1.403 


0.6 


1.430 


-3.2 


2.7 


2.1 


0.003 


1.405 


0.6 


1.422 


-3.2 


2.69 


2.7 



Table 2: The value of the amplitude of the jump of the temperature a obtained by analytical 
formula ([cr an ]) and by numerical simulation ([cx ca ;]). 



The functions rio(t), i = 1,2 are determined as the limits 

r i0 = ]im ri(t,e) 

in the asymptotic formulas. So in the capacity of 4fi"io{t*) we should take the constants which 
are obtained by differentiation of the functions ri(t,e) in the interval where this functions 
the most close to the line. 

We note that the positions of the free boundaries is fuzzy respect to e as e > (the width 
of the transition zone is proportional to e). Therefore to disclose the behavior of the shifting 
the free boundaries we consider the dynamic of the point of the crossing the order function 
u(r,t) and the axis r, i.e. we trace the dynamic of the points r(t) in which u(r(t),t) = for 
any fixed t, see FigflBl HU The instants of time t* in which the graphics of ri(t) and ^(t) 
become the horizontal lines correspond to the instants of time in which the order function 
u(r, t) becomes positive. This instants of time naturally do not equal to instants of time t min 
in which the amplitude of the jump of the temperature is maximal, see Table [TJ 

To calculate the amplitude of the jump of the temperature we determine the minimal 
value of the temperature a min corresponded to the value of the coordinate r min , see Table [H 
For example in Fig. [TUJ this value of the coordinate is become on the curve 4. The amplitude 
is equal to the absolute value of the difference between a m i n and the value of the temperature 
in the point r min and at the instant of time corresponded to the start of the fast changing 
of the temperature in this point (in Fig. [TU] the curve 3 gives this value of the temperature). 

The results of the verification of the formula (1181) are showed in Table [21 

In column [a an ] the data for the jump are given obtained by formula (j5J), and in column 
[(Teal] t ne data for the jump are given calculated by numerical simulation. From Table [2] we 
see that the analytical data converge to the numerical data as e decreases. 
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